logo

This script generates static and animated maps of prescribed burns in Three Rivers parks. It aligns the master burn spreadsheet with the spatial burn unit layer before plotting the maps. The maps can be used in annual reports and public presentations.

Setup

Load packages & data

# Loading required packages for analysis:
library(here)
library(tidyverse)
library(ggplot2)
library(sf)
library(ggspatial)
library(gganimate)
library(gifski)
library(transformr)
library(ggrepel)

source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
source(here("Scripts", "4_Covariate_Data_Prep.R"))
# source(here("Scripts", "5_Prep_Species.R"))

Data manipulation

The burn units spatial data are imported with the script “2_Spatial_Data_Import.R”. The burn data are imported and processed in the script “4_Covariate_Data_Prep.R”. We make copies of both of those layers for further modification here. There is lots of manual work to align both.

# Load spatial data on burn units
units <- spatial.mgmt.burns

# Load burn dates
burn_dates <- burns %>%
  select(Burn_unit_ID_new_2015, Burn_unit_description, year, burn_status_simple) %>%
  mutate(burn_status_simple = ifelse(is.na(burn_status_simple), 0, burn_status_simple))

# Manually make burn unit IDs in burn_dates match unit IDs in units
burn_dates <- burn_dates %>%
  mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "East Prairie", "NHTC 1", Burn_unit_ID_new_2015)) %>%
  mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "West Prairie", "NHTC 2", Burn_unit_ID_new_2015)) %>%
  mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_ID_new_2015 == "REB 1", "REB 1A", Burn_unit_ID_new_2015)) %>%
  mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_ID_new_2015 == "REB 2", "REB 1B", Burn_unit_ID_new_2015)) %>%
  mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "Open Area", "SPR NA", Burn_unit_ID_new_2015)) %>%
  # Carver 22 has west and east (a & b) in the burn dates layer, but just one polygon in the spatial layer. CAR 22 B in the burn dates layer has no burns marked. so I just use the A burns and apply them to the whole 22 polygon (both halves).
  filter(Burn_unit_ID_new_2015 != "CAR 22 B") %>% # drop 22 B
  mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "Old field south of mitigation pond. West Half", "CAR 22", Burn_unit_ID_new_2015)) %>%
  # Murphy 7 now has north and south in the burn dates layer (though both labeled 7B...), but just one polygon in the spatial layer. Going to combine both rows in the burns data and call both 7.
  filter(Burn_unit_description != "1996 Metro Prairie-East - sw marsh. South Half") %>% # drop South half
  mutate(burn_status_simple = ifelse(Burn_unit_ID_new_2015 == "MUH 7 B" & year == "2023", "1", burn_status_simple)) %>% # rescue the one burn in the south and apply to unit 7. TK this will need to be updated in future years (until spatial layer is updated...)
  mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_ID_new_2015 == "MUH 7 B", "MUH 7", Burn_unit_ID_new_2015))


# Join unit data (burn unit polygons) and burn data (data frame indicating when each unit was burned)
units <- units %>%
  select(ABBREVIATI, BURN_ID_2015, NRM_UNIT, Unit_ID) %>%
  rename(
    park = ABBREVIATI,
    id = BURN_ID_2015
  ) %>%
  # Alter some unit names/numbers to resolve differences between burn units & burn dates
  mutate(id = ifelse(park == "ELM_A", "3A", id)) %>%
  mutate(id = ifelse(park == "ELM_B", "3B", id)) %>%
  mutate(id = ifelse(park == "ELM_C", "3C", id)) %>%
  mutate(park = ifelse(park %in% c("ELM_A", "ELM_B", "ELM_C"), "ELM", park)) %>%
  mutate(id = ifelse(park == "GW" & NRM_UNIT == "GW1", "1", id)) %>%
  mutate(id = ifelse(park == "GW" & NRM_UNIT == "GW2", "2", id)) %>%
  mutate(id = ifelse(park == "GW" & NRM_UNIT == "GW3", "3", id)) %>%
  mutate(id = ifelse(park == "HYL  A" & id == "4", "4A", id)) %>%
  mutate(id = ifelse(park == "HYL  B" & id == "4", "4B", id)) %>%
  mutate(id = ifelse(park == "HYL  A" & id == "7", "7A", id)) %>%
  mutate(id = ifelse(park == "HYL  B" & id == "7", "7B", id)) %>%
  mutate(id = ifelse(park == "HYL  C" & id == "7", "7C", id)) %>%
  mutate(id = ifelse(park == "HYL" & NRM_UNIT == "B" & id == "6", "6A", id)) %>% # mis-match in NRM unit letter and renamed value intentional, pending confirmation [TK]
  mutate(id = ifelse(park == "HYL" & NRM_UNIT == "A" & id == "6", "6B", id)) %>% # mis-match in NRM unit letter and renamed value intentional, pending confirmation [TK]
  mutate(id = ifelse(park == "HYL  C" & id == "6", "6C", id)) %>%
  mutate(park = ifelse(park %in% c("HYL  A", "HYL  B", "HYL  C"), "HYL", park)) %>%
  mutate(park = ifelse(park %in% c("MUL_E", "MUL_W"), "MUL", park)) %>%
  mutate(id = ifelse(Unit_ID %in% c("NHTC1"), "1", id)) %>%
  mutate(id = ifelse(Unit_ID %in% c("NHTC2"), "2", id)) %>%
  mutate(id = ifelse(park == "REB A", "1A", id)) %>%
  mutate(id = ifelse(park == "REB B", "1B", id)) %>%
  mutate(park = ifelse(park %in% c("REB A", "REB B"), "REB", park)) %>%
  mutate(unit = paste0(park, " ", id)) %>% # Make match how units listed in burn data
  select(park, unit)

# Match units in spatial layer with names in date layer
units <- units %>%
  left_join(burn_dates, by = c("unit" = "Burn_unit_ID_new_2015")) %>%
  drop_na() %>%
  arrange(unit, year)

Generate static maps

Prep basemaps

We use habitat types as a basemap over which to plot the burns.

# Generate basemap (habitat types)
hab_groups_forest <- c("Oak Woods", "Early Successional Forest", "Shrubland", "Maple Basswood", "Flood Plain Forest", "Tamarack Bog")
hab_groups_grasslands <- c("Oldfield", "Prairie")
hab_groups_wetlands <- c("Wetland Cattail", "Wetland Sedge", "Wetland Lily", "Beach Mudflat")
hab_group_water <- c("Open Water")
hab_groups_cultivated <- c("Planted Cultivated")

basemap_habs <- spatial.habitats.current %>%
  mutate(type_simplified = ifelse(MLCCS_Data %in% hab_groups_forest, "forest",
    ifelse(MLCCS_Data %in% hab_groups_grasslands, "grassland",
      ifelse(MLCCS_Data %in% hab_groups_wetlands, "wetland",
        ifelse(MLCCS_Data %in% hab_groups_cultivated, "cultivated",
          ifelse(MLCCS_Data %in% hab_group_water, "water", NA)
        )
      )
    )
  )) %>%
  filter(!is.na(type_simplified))

basemap_cols <- c("forest" = "darkseagreen", "grassland" = "lightgoldenrod", "cultivated" = "wheat3", "wetland" = "lightblue", "water" = "skyblue4")

basemap <- ggplot() +
  geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE) +
  scale_fill_manual(values = basemap_cols) +
  theme_void()

basemap

Function for mapping individual parks/years

# Generate park name lookup table for parks with burns (relates codes in burn data to park names in park park boundaries dataset).
park_name_lookup <- data.frame(
  park = c("BAK", "BRY", "CAR", "CRO", "ELM", "FR", "GW", "HYL", "KW", "MUH", "MUL", "NHTC", "REB", "SPR"),
  Name = c("Baker", "Bryant", "Carver", "Crow-Hassan", "Elm Creek", "French", "Gale Woods", "Hyland", "Kingswood", "Murphy-Hanrehan", "The Landing", "NHTC", "Lake Rebecca", "Spring Lake Park")
)

# Function for quickly defining map extent in ggplot based on a particular park. Park name must match that listed in park boundaries sf.
map_extent_park <- function(park) {
  row <- which(spatial.boundaries.parks$Name == park)
  bbox <- st_bbox(spatial.boundaries.parks[row, ])
  coord_sf(
    xlim = c(bbox$xmin, bbox$xmax),
    ylim = c(bbox$ymin, bbox$ymax)
  )
}

# Function for mapping burns in a particular park/year
map_burns <- function(park, display_year) {
  units_filter_year <- filter(units, year == display_year)
  cols <- c("1" = "orangered3", "forest" = "darkseagreen", "grassland" = "lightgoldenrod", "wetland" = "lightblue", "water" = "skyblue4", cultivated = "wheat3")
  ggplot() +
    geom_sf(data = spatial.boundaries.parks, fill = "grey97", size = 0) +
    geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE, alpha = 0.7) +
    # geom_sf(data = spatial.infrastructure.trails %>%
    #               filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
    #               st_zm(drop = TRUE),
    #        linetype = "dotted", size = .5, alpha = .8, color = "tan4") +
    geom_sf(
      data = spatial.infrastructure.roads %>%
        filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
        st_zm(drop = TRUE),
      linetype = "solid", size = .4, alpha = .8, color = "tan4"
    ) +
    geom_sf(data = units_filter_year %>% filter(burn_status_simple == 1), aes(fill = as.factor(burn_status_simple)), show.legend = FALSE, color = "white", alpha = 0.75) +
    geom_sf(data = spatial.boundaries.parks, fill = NA, size = 0.7) +
    scale_fill_manual(values = cols) +
    scale_y_continuous(breaks = seq(40, 50, by = 0.02)) +
    scale_x_continuous(breaks = seq(-95, -90, by = 0.02)) +
    theme_minimal() +
    theme(
      axis.title = element_blank(),
      plot.title = element_text(size = 11, face = "bold", color = "orangered3"),
      plot.subtitle = element_text(size = 10, face = "plain")
    ) +
    map_extent_park(park) +
    annotation_scale(location = "br", style = "bar", bar_cols = c("grey40", "grey40"), line_col = "grey40", text_col = "grey40", line_width = 0.5) +
    labs(
      title = "Prescribed burns",
      subtitle = paste0(park, ": ", display_year)
    )
}

Example map, Crow-Hassan 1990:

map_burns("Crow-Hassan", "1990")

2023 burn maps

maps_2023 <- lapply(park_name_lookup$Name, map_burns, "2023")
names(maps_2023) <- park_name_lookup$Name
maps_2023
## $Baker

## 
## $Bryant

## 
## $Carver

## 
## $`Crow-Hassan`

## 
## $`Elm Creek`

## 
## $French

## 
## $`Gale Woods`

## 
## $Hyland

## 
## $Kingswood

## 
## $`Murphy-Hanrehan`

## 
## $`The Landing`

## 
## $NHTC

## 
## $`Lake Rebecca`

## 
## $`Spring Lake Park`

Generate animated maps

Function for animated maps

# Function for generating an animation of burns in a particular park
map_burns_for_ani <- function(park) {
  # Define colors
  cols <- c("1" = "orangered3", "cultivated" = "wheat3", "forest" = "darkseagreen", "grassland" = "lightgoldenrod", "wetland" = "lightblue", "water" = "skyblue4")

  # Pull burned untis only
  units_burned <- units %>% filter(burn_status_simple == 1)

  # Get park code from park name
  park_code <- park_name_lookup$park[which(park_name_lookup$Name == park)]

  # Calculate the minimum start year for the current park
  min_start_year <- min(as.numeric(as.character(units_burned$year[which(units_burned$park == park_code)])))

  # Create a subset of data starting from the year before the first burn
  units_subset <- units_burned %>% filter(park == park_code, as.numeric(as.character(year)) >= min_start_year - 1)

  # For mapping purposes, we need a fake/placeholder row in the birn units database (one for each unique combination of park and year). I'll have the geometry be a small circle centered on Minneapolis (out of frame of any park).
  fake_geometry <- st_sfc(st_point(c(-93.2650, 44.9778)))
  fake_geometry <- st_buffer(fake_geometry, dist = 100)
  st_crs(fake_geometry) <- st_crs(units)

  # Create a data frame with unique combinations of year and park
  fake_data <- expand.grid(year = unique(units$year), park = unique(units$park), burn_status_simple = 1, unit = NA, Burn_unit_description = NA) %>%
    filter(as.numeric(as.character(year)) >= min_start_year - 1) %>%
    filter(park == park_code)

  # Combine the data frame with the spatial data
  fake_df <- st_sf(fake_data, Shape = fake_geometry)

  # Bind the fake rows to the original units_burned dataframe
  units_subset <- rbind(units_subset, fake_df)

  # Generate basemap for animation
  plot <- ggplot() +
    geom_sf(data = spatial.boundaries.parks, fill = "grey97", size = 0) +
    geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE, alpha = 0.7) +
    # geom_sf(data = spatial.infrastructure.trails %>%
    #               filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
    #               st_zm(drop = TRUE),
    #        linetype = "dotted", size = .5, alpha = .8, color = "tan4") +
    geom_sf(
      data = spatial.infrastructure.roads %>%
        filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
        st_zm(drop = TRUE),
      linetype = "solid", size = .4, alpha = .8, color = "tan4"
    ) +
    geom_sf(data = units_subset, aes(fill = as.factor(burn_status_simple)), show.legend = FALSE, color = "white", alpha = 0.75) +
    geom_sf(data = spatial.boundaries.parks, fill = NA, size = 0.7) +
    scale_fill_manual(values = cols) +
    scale_y_continuous(breaks = seq(40, 50, by = 0.02)) +
    scale_x_continuous(breaks = seq(-95, -90, by = 0.02)) +
    theme_minimal() +
    theme(
      axis.title = element_blank(),
      plot.title = element_text(size = 11, face = "bold", color = "orangered3"),
      plot.subtitle = element_text(size = 10, face = "plain")
    ) +
    map_extent_park(park) +
    annotation_scale(location = "br", style = "bar", bar_cols = c("grey40", "grey40"), line_col = "grey40", text_col = "grey40", line_width = 0.5) +
    labs(
      title = "Prescribed burns",
      subtitle = paste0(park, ": {next_state}")
    ) # not sure why next_state works but closest_State doesn't; closest lags the data by 1 year

  # Facet by year for animation frames
  plot_animation <- plot + transition_states(
    year,
    transition_length = 0,
    state_length = 1,
    wrap = TRUE
  )

  # Generate animation
  mapGIF <- animate(
    plot_animation,
    nframes = length(unique(units_subset$year)) + 2,
    fps = 1.5,
    end_pause = 2,
    renderer = gifski_renderer()
  )

  # Save animation
  anim_save(
    filename = paste0("burns_", park, ".gif"),
    animation = mapGIF,
    path = here("Results", "990_burn_animation")
  )
}

Generate & save animations

parks_list <- c("Baker", "Bryant", "Carver", "Crow-Hassan", "Elm Creek", "French", "Gale Woods", "Hyland", "Kingswood", "Murphy-Hanrehan", "The Landing", "Lake Rebecca", "Spring Lake Park") # parks with  burn data
lapply(parks_list, map_burns_for_ani)

View animations

files <- list.files(path = here("Results", "990_burn_animation"), pattern = ".gif", full.names = TRUE)
files_names <- list.files(path = here("Results", "990_burn_animation"), pattern = ".gif", full.names = FALSE)
files_names <- gsub("\\.gif", "", files_names)
files_names <- gsub("burns_", "", files_names)

for (i in 1:length(files)) {
  cat(paste0("#### ", files_names[i], "\n"))
  cat(paste0("![](", files[i], ")\n\n"))
}

Baker

Bryant

Carver

Crow-Hassan

Elm Creek

French

Gale Woods

Hyland

Kingswood

Lake Rebecca

Murphy-Hanrehan

Spring Lake Park

The Landing

 




By Sam Safran